namespace NBody.Analysis
{
    using System;
    using Gsl;
    
    public class Shape
    {
        public static void CalculateGlobal(ref NBodySystem sys, out double q, out double s, double tolerance)
        {   
            Console.Error.WriteLine("Calculating global shape of system");
            double old_q, old_s;         
            q = s = 1.0;
            
            // Iterative procedure.  See Dubinski and Carlberg (1991).
            Matrix m = new Matrix(3, 3);
            Console.Error.Write("Error = {0:e3}", 0.0);
            Gsl.Eigensystems.RealSymmetric eigensystem = new Gsl.Eigensystems.RealSymmetric(3);
            do
            {
                m.SetZero();
                // calculate the modified inertial tensor
                for (int i = 0; i < sys.NumParts; i++)
                {
                    double a2 = sys[i].X * sys[i].X + sys[i].Y * sys[i].Y / q / q + sys[i].Z * sys[i].Z / s / s;
                    for (int j = 0; j < 3; j++)
                    {
                        for (int k = 0; k < 3; k++)
                        {
                            m[j, k] += sys[i].Position(j) * sys[i].Position(k) / a2;
                        }
                    }
                }
                
                Vector values;
                Matrix vectors;
                eigensystem.Calculate(ref m, out values, out vectors);
            
                old_q = q;
                old_s = s;
                        
                q = Math.Sqrt(values[1] / values[0]);
                s = Math.Sqrt(values[2] / values[0]);
            
                RotateParticles(ref sys, vectors);
            
                Console.Error.Write("\b\b\b\b\b\b\b\b\b\b{0:e3}", Math.Abs(old_s - s));
            } while (Math.Abs(old_s - s) > tolerance || Math.Abs(old_q - q) > tolerance);
            Console.Error.WriteLine();            
        }
        
        private static void RotateParticles(ref NBodySystem sys, Matrix vectors)
        {
            double x, y, z;
            double vx, vy, vz;
            for (int i = 0; i < sys.NumParts; i++)
            {
                x = y = z = 0.0;
                vx = vy = vz = 0.0;
                for (int j = 0; j < 3; j++)
                {
                    x += vectors[0, j] * sys[i].Position(j);
                    y += vectors[1, j] * sys[i].Position(j);
                    z += vectors[2, j] * sys[i].Position(j);
                
                    vx += vectors[0, j] * sys[i].Velocity(j);
                    vy += vectors[1, j] * sys[i].Velocity(j);
                    vz += vectors[2, j] * sys[i].Velocity(j);
                }
                sys[i].SetPosition((double)x, (double)y, (double)z);
                sys[i].SetVelocity((double)vx, (double)vy, (double)vz);
            }
        }
    }
}
            
